st_layers(here("prac3_data", "gadm36_AUS.gpkg"))
Driver: GPKG
Available layers:
Ausoutline <- st_read(here("prac3_data", "gadm36_AUS.gpkg"),
layer='gadm36_AUS_0')
Reading layer `gadm36_AUS_0' from data source `D:\CASA\GIS\prac3_data\gadm36_AUS.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS: WGS 84
# re-project or transform CRS
AusoutlinePROJECTED <- Ausoutline %>%
st_transform(.,3112) # GDA94, a local CRS for Australia
print(AusoutlinePROJECTED)
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -2083066 ymin: -6460625 xmax: 2346599 ymax: -1115948
Projected CRS: GDA94 / Geoscience Australia Lambert
GID_0 NAME_0 geom
1 AUS Australia MULTIPOLYGON (((1775780 -64...
##From sf to sp
#AusoutlineSP <- Ausoutline %>%
# as(., "Spatial")
##From sp to sf
#AusoutlineSF <- AusoutlineSP %>%
# st_as_sf()
library(raster)
library(terra)
jan<-terra::rast(here("prac3_data", "wc2.1_5m_tavg_01.tif"))
# have a look at the raster layer jan
jan
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : wc2.1_5m_tavg_01.tif
name : wc2.1_5m_tavg_01
min value : -46.042
max value : 34.065
plot(jan)
using the Mollweide projection saved to a new object. The Mollweide projection retains area proportions whilst compromising accuracy of angle and shape
# set the proj 4 to a new object
pr1 <- terra::project(jan, "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
#or....
newproj<-"ESRI:54009"
# get the jan raster and give it the new proj4
pr1 <- jan %>%
terra::project(., newproj)
plot(pr1)
#back to WGS 84
pr1 <- pr1 %>%
terra::project(., "EPSG:4326")
plot(pr1)
# look in our folder, find the files that end with .tif and
library(fs)
dir_info("D:/CASA/GIS/prac3_data/")
# select the data we actually want
library(tidyverse)
listfiles<-dir_info("D:/CASA/GIS/prac3_data/") %>%
filter(str_detect(path, ".tif")) %>%
dplyr::select(path)%>%
pull()
#have a look at the file names
listfiles
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_01.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_02.tif
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_03.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_04.tif
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_05.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_06.tif
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_07.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_08.tif
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_09.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_10.tif
D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_11.tif D:/CASA/GIS/prac3_data/wc2.1_5m_tavg_12.tif
Then load all of the data straight into a SpatRaster. A SpatRaster is a collection of raster layers with the same spatial extent and resolution.
worldclimtemp <- listfiles %>%
terra::rast()
#have a look at the raster stack
worldclimtemp
class : SpatRaster
dimensions : 2160, 4320, 12 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
sources : wc2.1_5m_tavg_01.tif
wc2.1_5m_tavg_02.tif
wc2.1_5m_tavg_03.tif
... and 9 more source(s)
names : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ...
min values : -46.042, -44.800, -57.986, -64.200, -64.829, -64.395, ...
max values : 34.065, 32.908, 33.081, 34.277, 36.299, 38.458, ...
# access the january layer
worldclimtemp[[1]]
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : wc2.1_5m_tavg_01.tif
name : wc2.1_5m_tavg_01
min value : -46.042
max value : 34.065
#rename the layers
month <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
names(worldclimtemp) <- month
# get data for January using new layer name
worldclimtemp$Jan
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : wc2.1_5m_tavg_01.tif
name : Jan
min value : -46.042
max value : 34.065
Using a raster stack we can extract data with a single command!! For example, make a dataframe of some sample sites — Australian cities/towns.
site <- c("Brisbane", "Melbourne", "Perth", "Sydney", "Broome", "Darwin", "Orange",
"Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle",
"Wollongong", "Logan City" )
lon <- c(153.03, 144.96, 115.86, 151.21, 122.23, 130.84, 149.10, 115.64, 145.77,
138.6, 153.43, 149.13, 151.78, 150.89, 153.12)
lat <- c(-27.47, -37.91, -31.95, -33.87, 17.96, -12.46, -33.28, -33.33, -16.92,
-34.93, -28, -35.28, -32.93, -34.42, -27.64)
#Put all of this inforamtion into one list
samples <- data.frame(site, lon, lat, row.names="site")
# Extract the data from the Rasterstack for all points
AUcitytemp<- terra::extract(worldclimtemp, samples)
Add the city names to the rows of AUcitytemp
Aucitytemp2 <- AUcitytemp %>%
as_tibble()%>%
add_column(Site = site, .before = "Jan")
take Perth as an example. We can subset our data either using the row name:
Perthtemp <- Aucitytemp2 %>%
filter(site=="Perth")
Make a histogram of Perth’s temperature. The tibble stored the data as double and the base hist() function needs it as numeric..
hist(as.numeric(Perthtemp))
Warning: NAs introduced by coercion
library(tidyverse)
#define where you want the breaks in the historgram
userbreak<-c(8,10,12,14,16,18,20,22,24,26)
# remove the ID and site columns
Perthtemp <- Aucitytemp2 %>%
filter(site=="Perth")
t<-Perthtemp %>%
dplyr::select(Jan:Dec)
hist((as.numeric(t)),
breaks=userbreak,
col="red",
main="Histogram of Perth Temperature",
xlab="Temperature",
ylab="Frequency")
#Check out the histogram information R generated
histinfo <- as.numeric(t) %>%
as.numeric()%>%
hist(.)
histinfo
$breaks
[1] 12 14 16 18 20 22 24 26
$counts
[1] 2 2 2 1 1 2 2
$density
[1] 0.08333333 0.08333333 0.08333333 0.04166667 0.04166667 0.08333333 0.08333333
$mids
[1] 13 15 17 19 21 23 25
$xname
[1] "."
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
breaks — the cut off points for the bins (or bars), we just specified these counts — the number of cells in each bin midpoints — the middle value for each bin density — the density of data per bin
Check the layer by plotting the geometr
plot(Ausoutline$geom)
#simplify the .shp with lots of points
AusoutSIMPLE <- Ausoutline %>%
st_simplify(., dTolerance = 1000) %>% #controls the level of generalisation in the units of the map
st_geometry()%>%
plot()
make sure that both of our layers are in the same coordinate reference system before combine
print(Ausoutline)
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS: WGS 84
GID_0 NAME_0 geom
1 AUS Australia MULTIPOLYGON (((158.6928 -5...
#this works nicely for rasters
crs(worldclimtemp)
[1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"
Austemp <- Ausoutline %>%
# now crop our temp data to the extent
terra::crop(worldclimtemp,.)
# plot the output
plot(Austemp)
exactAus<-terra::mask(Austemp, Ausoutline)
exactAus
class : SpatRaster
dimensions : 551, 554, 12 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : 112.9167, 159.0833, -55.08333, -9.166667 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : memory
names : Jan, Feb, Mar, Apr, May, Jun, ...
min values : 6.4875, 6.53125, 5.845714, 4.571428, 2.642, -0.394, ...
max values : 34.0650, 32.90800, 31.986000, 30.400000, 28.800, 27.200, ...
#subset using the known location of the raster
hist(exactAus[[3]], col="red", main ="March temperature")
#make our raster into a data.frame to be compatible with ggplot2, using a dataframe or tibble
exactAusdf <- exactAus %>%
as.data.frame()
library(ggplot2)
# set up the basic histogram
gghist <- ggplot(exactAusdf,
aes(x=Mar)) +
geom_histogram(color="black",
fill="white")+
labs(title="Ggplot2 histogram of Australian March temperatures",
x="Temperature",
y="Frequency")
# add a vertical line to the hisogram showing mean tempearture
gghist + geom_vline(aes(xintercept=mean(Mar,
na.rm=TRUE)),
color="blue",
linetype="dashed",
size=1)+
theme(plot.title = element_text(hjust = 0.5))
put our variable (months) into a one column using pivot_longer()
squishdata<-exactAusdf%>%
pivot_longer(
cols = 1:12,
names_to = "Month",
values_to = "Temp"
)
select two month
twomonths <- squishdata %>%
# | = OR
filter(., Month=="Jan" | Month=="Jun")
meantwomonths <- twomonths %>%
group_by(Month) %>%
summarise(mean=mean(Temp, na.rm=TRUE))
meantwomonths
ggplot(twomonths, aes(x=Temp, color=Month, fill=Month)) +
geom_histogram(position="identity", alpha=0.5)+
geom_vline(data=meantwomonths,
aes(xintercept=mean,
color=Month),
linetype="dashed")+
labs(title="Ggplot2 histogram of Australian Jan and Jun
temperatures",
x="Temperature",
y="Frequency")+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
data_complete_cases <- squishdata %>%
drop_na()%>% #dropped all the NAs
# Month column map in descending order (e.g. Jan, Feb, March..)
mutate(Month = factor(Month, levels = c("Jan","Feb","Mar",
"Apr","May","Jun",
"Jul","Aug","Sep",
"Oct","Nov","Dec")))
# Plot faceted histogram
ggplot(data_complete_cases, aes(x=Temp, na.rm=TRUE))+
geom_histogram(color="black", binwidth = 5)+
labs(title="Ggplot2 faceted histogram of Australian temperatures",
x="Temperature",
y="Frequency")+
facet_grid(Month ~ .)+
theme(plot.title = element_text(hjust = 0.5))
# mean per month
meanofall <- squishdata %>%
group_by(Month) %>%
summarise(mean = mean(Temp, na.rm=TRUE))
# print the top 1
head(meanofall, n=1)
# standard deviation per month
sdofall <- squishdata %>%
group_by(Month) %>%
summarize(sd = sd(Temp, na.rm=TRUE))
# maximum per month
maxofall <- squishdata %>%
group_by(Month) %>%
summarize(max = max(Temp, na.rm=TRUE))
# minimum per month
minofall <- squishdata %>%
group_by(Month) %>%
summarize(min = min(Temp, na.rm=TRUE))
# Interquartlie range per month
IQRofall <- squishdata %>%
group_by(Month) %>%
summarize(IQR = IQR(Temp, na.rm=TRUE))
# perhaps you want to store multiple outputs in one list..
lotsofstats <- squishdata %>%
group_by(Month) %>%
summarize(IQR = IQR(Temp, na.rm=TRUE),
max=max(Temp, na.rm=T))
# or you want to know the mean (or some other stat)
#for the whole year as opposed to each month...
meanwholeyear=squishdata %>%
summarize(meanyear = mean(Temp, na.rm=TRUE))
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.